Plug flow and the breakdown of Bagnold scaling in cohesive granular flows 
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Cohesive granular media flowing down an inclined plane are studied by discrete element simu- 
lations. Previous work on cohesionless granular media demonstrated that within the steady flow 
regime where gravitational energy is balanced by dissipation arising from intergrain forces, the ve- 
locity profile in the flow direction scales with depth in a manner consistent with the predictions of 
Bagnold. Here we demonstrate that this Bagnold scaling does not hold for the analogous steady- 
flows in cohesive granular media. We develop a generalization of the Bagnold constitutive relation 
to account for our observation and speculate as to the underlying physical mechanisms responsible 
for the different constitutive laws for cohesive and noncohesive granular media. 

PACS numbers: 45.70.-n, 45.70.Mg, 83.10.Ff 



I. INTRODUCTION 

One of the central questions in the study of granu- 
lar flows is how to determine the relationship between 
the microphysics of grain interactions and the collective 
or macroscopic flow properties of the system. Elucidat- 
ing these flow properties is of fundamental importance 
to a variety of fields ranging from civil engineering to 
geophysics Moreover, this question remains at the 
forefront of many-body physics as its solution appears 
to demand entirely new concepts that are applicable to 
systems driven far from equilibrium. 

In this article we explore the velocity field in steadily 
flowing granular media in the chute geometry - an in- 
clined plane having a rough surface at its base and a free 
surface at the top. Such flows have beenpreviously char- 
acterized both experimentally 1111 I I |H H S 13 and 
via simulation jTfl, EH El EI Oil El In this work 
we address, through large-scale discrete element simula- 
tions, the change in the flow field in the material as a 
function of interparticle adhesion. We present a modi- 
fied constitutive law relating shear stress within the bulk 
to the local gradients in the velocity field and discuss 
the implications of our proposed continuum description 
of the material. In addition we observe the formation of 
a plug flow regime in the cohesive slab extending from 
the free surface at the top of the slab into the bulk. The 
depth of this plug flow regime can be calculated from a 
comparison of the yield stress of the cohesive material 
to depth-dependent shear stress due to the weight of the 
slab. 

The study of cohesive granular media is necessary for 
the application of granular physics to the behavior of such 
materials in their geophysical context, i.e. the dynamics 
of landslides, snow avalanches, and even the lunar re- 
golith. In the former examples, a small aqueous wetting 
layer on the particles forms menisci at intergrain con- 
tacts to produce an adhesive force. In the latter case, 
high-vacuum conditions facilitate close contacts between 



particles leading to interparticle adhesion due to van der 
Waals interactions. The study of cohesive granular ma- 
terials also helps to elucidate fundamental questions con- 
cerning the physical description of this state of matter. 
One aspect of granular materials that differentiates them 
from random elastic solids or liquids is that cohesionless 
granular system cannot support tensile stresses. By ob- 
serving the quantitative effect of known cohesive forces 
upon granular flows, one should be able to bridge be- 
tween the better understood elastic solids and/or viscous 
liquids and granular materials, which, by comparison, re- 
main rather mysterious. The change in the stress-strain 
relation due to internal cohesion also provides insight into 
the rheology of cohesionless granular matter. We discuss 
this further below. 

In the study of granular mechanics it is clearly de- 
sirable to develop a continuum description of flow since 
calculating the detailed dynamics of numerous inter- 
grain collisions rapidly becomes intractable with increas- 
ing numbers of such particles. In addition one expects, 
based on experience with the continuum mechanics of 
solids, liquids, and gases, to be able to develop a set of 
relations between macroscopic averaged observable such 
as the velocity field u Q (x), mass density p(x), and stress 
tensor a a p (x) within the steadily- flowing granular mate- 
rial where details of the grai n interactions enter through 
a small set of parameters [rH Hil IT^ . Such a constitutive 
relation between the stress state in the material and its 
rate of deformation taken together with momentum con- 
servation completely determines the flow properties of 
the granular system assuming the boundary conditions 
on the flow are sufficiently well-known. There are now 
a number of such hydrodynamic descriptions of granu- 
lar flow which seek to derive such constitutive equations 
from a more microscopic theory [l7L l20| . There are a 
variety of such models that rely on arguments based on 
effective viscosities 0,111 |23j, transient force-chains p4) 
or a superposition of a rate-dependent contribution aris- 
ing from collisional interactions and a rate-independent 
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part related to enduring frictional contacts among the 
grains p5|. 

A well-known proposal by Bagnoldj^ for such a con- 
stitutive relation for granular flows is that the shear stress 
in the flow is proportional to the instantaneous square of 
the rate of strain tensor and that the density is constant 
throughout the material. Taking (as we do throughout 
this work) a coordinate system where the a;-axis is di- 
rected in the flow direction and the z-axis as the upward 
normal to the inclined plane, the Bagnold relation is 

<Jxz = k (d z v x ) 2 , (1) 

where n is a constant having the dimensions of M/L. 
One can reduce the supporting argument for Eq. ^ to 
dimensional analysis; in the flowing granular system the 
one energy density available is the kinetic energy of the 
grains themselves ~ pv 2 . As long as the local mass den- 
sity remains constant, Galilean invariance requires the 
shear stress to be proportional to the square of the veloc- 
ity difference across the sample. The further assumption 
that the relationship between stress and the rate of strain 
is local produces Eq. ^ A simple heuristic argument to 
derive the scaling relation shown in Eq.^and attributed 
to Bagnold is as follows: the transport in the gradient 
direction (z) of the component of momentum parallel to 
the velocity direction (x) occurs only through collisions 
between grains. The rate of those collisions depends on 
the average velocity difference between grains at different 
locations along the gradient direction ~ ad z v x where a 
is the grain radius. The momentum transfer per collision 
also scales linearly with the velocity difference leads to 
Cxz ~ m(ad z v x ) where m is the mass of a single grain. 

It is instructive to compare the above argument to the 
standard expression for the viscosity of a gas as deter- 
mined by kinetic theory. In that case the rate of inter- 
molecular collisions is controlled by the root mean square 
velocity of the gas molecules that is fixed by the equipar- 
tition theorem. Thus the viscosity of the gas is propor- 
tional to the square root of temperature. The granu- 
lar system is effectively at zero temperature so that the 
random particle velocities are driven themselves by the 
macroscopic velocity scale making the effective granular 
viscosity as extracted from Eq. ^ depend on shear rate. 

In the remainder of this article we first discuss the sim- 
ulation method in sectionllllbefore turning to a discussion 
of our numerical results and calculations in section II I II 
In that section we discuss the effect of intergrain adhe- 
sion on the validity of the Bagnold constitutive law. We 
then conclude in section llVl 

II. SIMULATION METHOD 

We performed discrete element simulations on a three- 
dimensional system of N monodisperse particles of mass 
m and diameter d on a rough base tilted an angle 9 with 
respect to gravity. Our simulations volume is a rectan- 
gular box with periodic boundary conditions in both the 




FIG. 1: (color online) A representative example of the flowing 
granular slab of height H — 300. In this cohesive granular slab 
(A = 0.8) the system separates into a solid-like plug (blue) 
sliding upon a flowing granular bed (red). The rough surface 
of inclined plane is shown in green. 



x and y directions, a rough base and an open top. The 
rough base is created by taking a slice through a pre- 
viously random close packed state of particles with the 
same diameter d. We define the z-axis to be normal to 
the base and the ir-axis as the direction of flow. We study 
a system of length 40c? and 10c? in the x and y directions 
respectively. We studied three granular slabs of differing 
heights H, containing 42000, 83000, and 125000 parti- 
cles respectively. We refer to these three systems by their 
approximate height, namely H = 100c?, 200c?, and 300c?. 
We show an example of a cohesive granular material in 
the steadily flowing state in Fig. ^ for H — 300c?. 

We employ a modified version of the model developed 
by Cundall and Strack to model cohesionless partic- 
ulates. The model uses Hookean contacts to model grain 
to grain interactions. We add a cohesive force between 
particles to simulate damp granular media. Two contact- 
ing spheres i and j with positions and Yj are separated 
by tij = Tj — Tj, the compression is then 5ij = d — |ry|. 
The relative velocity is v^- = Vj — Vj which can be sepa- 
rated into normal and tangential components 

v n tJ = (vy • n, 3 )n i3 (2) 

V t y = Vij - V„y - + UJj) X T^. (3) 

The normal and tangential forces acting on particle i due 
to particle j can then be written as 

vn 

F nij = (knSijiHj - y7nV„ y )+P^ (Sij) (4) 

771 

F t y = (-fetUiy - — TtVty), (5) 
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where k n ,k t and 7„,7t are the elastic and viscoelastic 
constants for interparticle motion along (n) and normal 
to the line of centers (t). u tij is the elastic tangential 
displacement that is set equal to zero at the initiation 
of contact and is truncated to satisfy the Coulomb yield 
criterion, |Ft 4 J < /x|F Wi J. The additional normal force 
Fy (6ij), is the cohesive normal force between particles i 
and j. This normal force is derived from an effective co- 
hesive potential acting on particle i due to particle j. We 
simply choose a gaussian well centered around particle j 

U^ = -Ae^, (6) 

where A is the strength of the attraction (in units of 
mgd) and I is the width of the well. The force is then 
P|. = -Vf/P- We focus on 0.0 < A < 1.0, where A = 1.0 
corresponds to a force capable of supporting 85 particles 
under gravity. While this form of the cohesive poten- 
tial cannot be justified in terms of details of either van 
der Waals or wetting-layer mediated interparticle inter- 
actions, we choose this form so that we can vary the 
strength of the interparticle cohesion via the well depth 
A while maintaining the short-range nature of the in- 
teraction, which is controlled by I. The force captures 
the essential features of more physical adhesion scenar- 
ios and we expect that the effects of interaction cohesion 
as modelled here, to accurately reflect the behavior of 
experimentally accessible cohesive granular media. 

Most of the simulations were run with k n = 2 x 10 5 ^p, 
k t = |fc„, 7„ = SOr^ 1 , and jt = 0, where r = yj d/g. 
For the case of Hookean springs used here, these param- 
eters give a coefficent of restitution e n — 0.88 force nor- 
mal collisions. The normal spring constant k n is large 
enough to minimize the overlap between particles but 
small enough to be computationally efficient. Previous 
simulations for cohcsionless grains found that this 
value of k n provides results similiar to those for larger 
values of k n . However we do find a subtle dependence of 
the flow rheology upon the value of k n which we explore 
later in this article. For the coefficent of friction, we use 
fi = 0.5, which is typical for the types of materials we are 
modeling. The normal damping term operates only when 
the particles are in actual contact so that, in a narrow 
window of interparticle separations d < r < d + £ these 
particles will experience a cohesive force without the ve- 
locity dependent damping. Interpreted in terms of ac- 
tion of interparticle wetting layers, this suggests that we 
discount dissipative hydrodynamics in the wetting layer 
over microscopic plastic deformation in the particle them- 
selves. We expect that such hydrodynamic effects modify 
the effective value of j n as well as the interparticle sepa- 
ration over which it operates, but we have not attempted 
to model this particular cohesion scenario in detail. 

In a gravitational field g, the total force on a particle 
from Eqs.^El and0 

Ff t = m i g + ^(F n ..+F t ..+F^.), (7) 
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where the sum is over neighboring particles. 

The stress tensor within a volume V is computed by 
summing over both the contact and kinetic terms of each 
particle within that volume 

i i=£j 

where F^- = Ff. . + F@. . + F^ and v is the time-averaged 
velocity of the particles in V. 

The timestep for the integration of the equations of 
motion is St = 10~ 4 r. After equilibration the systems 
were typically run between 1 — 5 x 10 7 time steps. Steady 
state was determined using the total kinetic energy of the 
system as a criterion for suitable equilibration of each 
sample. 



III. RESULTS 

In chute flow there are three qualitatively distinct 
regimes determined by the height of the slab H, the in- 
clination angle 9 of the lower surface with respect to the 
direction of gravity 0, ^| , and the intergrain cohesive 
stresses determined in our model by A - see Eq. |SJ For 
small enough inclination angles (i.e. below the maximum 
critical angle) or alternatively for short and/or cohesive 
slabs, the granular heap is stationary. Below the angle of 
repose for a given slab height and cohesive energy, tran- 
siently flowing states of the slab dissipate energy more 
rapidly than the input of gravitational potential energy 
so that the slab stops. We refer to this as the no flow 
regime. At much larger angles, on the other hand, the 
energy dissipation with the slab is less than the gravi- 
tational energy input so that the slab continuously ac- 
celerates down the plane. We refer to the parameter 
space exhibiting this behavior as the unstable regime. 
At angles intermediate between these two regimes for a 
given cohesive energy and slab height, we observe steady- 
state flows. In this work we concentrate primarily on 
such steady-state behavior in slabs of the largest heigth 
H = 300c?. The phase diagram spanned by the cohesion 
parameter A and the inclination angle of the slab 9 for 
this system shown in Fig. In this work we investigate 
the steady-state flowing regime after transient behavior 
associated with the initiation of flow have decayed. The 
onset of such steady-state behavior is determined by ob- 
serving the total kinetic energy of the system. 



A. Plug-flow 

In Fig. [3^ w e plot the velocity in the flow direction vs. 
height from the bottom of the slab (z) for a range of val- 
ues of the cohesive energy A. In each case the inclination 
angle is 9 = 22°. One of the central results of this work is 
immediately evident in these figures. For nonzero values 
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FIG. 2: Phase diagram for chute flow of a fixed height gran- 
ular slab. The diagram is spanned by the tilt angle 9 and 
strength of cohesion A. There exist three well defined regions 
corresponding to no flow, stable flow and unstable flow. Lines 
are drawn to seperate the regions. 



of the intergrain cohesive energy A, the velocity plateaus 
at a finite fraction of the height of the slab. The granular 
material has, in effect, phase separated into a liquid-like 
flowing region and a solid-like region that does not admit 
a non-zero rate of the shear strain. We refer to this be- 
havior near the free surface as plug flow and the solid-like 
region as the plug. It is clear from this figure that the 
thickness of the plug depends on the cohesive energy A. 
Figure 03 shows the dependence on the plug thickness 
at fixed A on the inclination angle 0. Finally Fig. 
demonstrates that the thickness of the plug is indepen- 
dent of the total height of the slab as long as the material 
remains in the flowing state for that slab height. Therein 
the granular velocity down the chute for three slabs of 
heights H = 100rf, 200<i, and 300d is plotted as a func- 
tion of the depth from the free surface for an inclination 
angle of 6 = 22° and A = 0.8. 

The development of the plug in steady state is signalled 
not only by the appearance of the velocity plateau, but 
also by the development of the localized jump in the par- 
ticle volume fraction with height. The flowing material 
below due to shear dilatancy is generically at a lower 
volume fraction, while the plug is denser. Typically this 
volume fraction change is on the order of 4% and occurs 
over a height range of 4 — 10 particle diameters. We use 
this abrupt change in the particle volume fraction with 
height to more precisely determine the interface between 
the flowing and plug regimes. The inset of Fig. 21 shows 
a typical example of this volume fraction change at the 
plug boundary. 

To interpret these observations, we propose that the 
cohesive granular material can support a finite yield 
stress a c before flowing. Since shear flow requires di- 
latancy, the critical yield stress should be equal to the 
typical cohesive stress in the pile. We estimate the max- 
imal value of this cohesive stress as follows. We assume 
all contacts provide some average adhesive force and then 
calculate the mean number of such interparticle contacts 
per unit cross-sectional area, n c , in terms of the volume 
fraction of the particles, 0. The yield stress is then esti- 
mated as the product of the average adhesive force per 
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FIG. 3: Velocity profiles as a function of z. In all three figures 
the velocity is measured in units of vo — \fgd. In (a) profiles 
are shown for different values of A with 9 — 22°, H — 300d. 
In (b) velocity profiles for different values of the tilt angle 9 
with A = 0.6, H = 300d. And (c) the velocity profile as a 
function of the depth from the free surface for differing slab 
heights H with A = 0.6, 9 = 22°. For H = lOOd, only a 
stopped plug is observed as H is similar to the size of the 
plug, w. In each profile two distinct regions are visible, a 
flowing state at depth and a plug flow of varying thicknesses 
at the free surface. For the noncohesive case A = 0.0 the plug 
size vanishes. 



contact and n c . To determine the average adhesive force 
per contact, /, we assume that the interparticle sepa- 
rations at each contact are randomly distributed within 
the attractive potential well between d and d + £. Using 
this assumption we find that the mean magnitude of the 
adhesive force is given by 



/ 



\U(0)-U(£)\ 



(9) 



in terms of the potential U defined in Eq. H3 From Eq. El 
and n c — 3</>/(27ra 2 ) we find that the yield stress is 
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(10) 



We determine the thickness of the plug w by equating the 
shear stress in the slab due to the weight of the overlying 
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FIG. 4: The thickness of the plug w vs. the cohesive energy 
A for a fixed angle of inclination, 9 — 22°, and pile height, 
H=300d. The dashed line shows the plug width prediction 
from Eq. 1111 The uncertainties in the plug thickness are de- 
termined from the width of the density jump at the lower 
boundary of the plug. Note that the plug boundary becomes 
more diffuse for smaller values of A. 



material to this yield stress. We find the thickness of the 
plug is 
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where e is the base of the natural logarithm. To test 
this proposal we plot in Fig. 0] the calculated plug from 
Eg. 1111 with that measured by particle density change at 
the plug boundary. 

Over the parameter range tested, the calculated thick- 
ness of the plug is in reasonable agreement with the data. 
Clearly the linear dependence of the plug size on the co- 
hesive energy is supported by the data. Precise tests of 
the 1 / sin 9 dependence of the plug thickness, on the other 
hand, are difficult due to the limited range of available 
angles for which the slab exhibits steady-state flow. Our 
calculation of the slope of the theoretical curve must be 
treated as an upper bound as it is based on the assump- 
tion that all interparticle contacts generate some attrac- 
tive force. Due to jamming in disordered sphere pack- 
ings, it is reasonable to suppose that at least some con- 
tacts have interparticle spacings less than d and are there- 
fore repulsive in nature. Such contacts between jammed 
spheres reduce the effective yield stress. 



B. Testing Bagnold-scaling 

From the Bagnold conjecture, Eq. 2] for the constitu- 
tive relation in the flowing granular state, one may imme- 
diately determine the flow profile down the chute. The 
resulting velocity profile takes the form: 



v x (z) 



H 3 pg sin 6 
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H 



3/2' 



(12) 



In Appendix^we derive this velocity profile for the case 
of a free surface at the top of the flowing medium, i. e. no 
plug as well as the analogous velocity profile for the case 



of a finite shear stress applied to the top of the flowing 
state due to the weight of the plug. 

In previous research on cohcsionlcss granular chute 
flow (A — 0.0), the down-chute velocity versus height 
v x (z) has been reported as being consistent with the Bag- 
nold power-law form. We also find such apparently rea- 
sonable agreement - see the uppermost curve in Fig. 
In the presence of cohesion (A > 0), however, the ex- 
pected Bagnold velocity profile fits the data poorly. We 
interpret this failure of the Bagnold hypothesis as demon- 
strating a new mode for vertically transporting (z) down- 
chute momentum (p x ) due to the presence of long-lived 
contacts in the material. These long-lived adhesive con- 
tacts in the material create clusters spanning different 
streamlines in the flow as has been proposed by Erta§ and 
Halsey [H| ■ In the presence of shearing flow, these forces 
acting along these clusters of particles transmit momen- 
tum p x proportional to the shear rate 7 = dv x /dz. Thus, 
for cohesive granular materials the Bagnold relation can 
be generalized to a form that is a sum of terms that are 
linear and quadratic in the shear rate. We propose the 
modified Bagnold relation, 

<J X z - <J C = k(— — ) +/3( — ). (13) 
oz oz 

The second of these terms represents the new mode of 
momentum transport made possible through the long- 
lived contact networks in the material while the first term 
arises from the short time scale collisions originally con- 
sidered by Bagnold [27]]. The constant stress term on 
the LHS of the above equation is the finite yield stress 
of the cohesive material. The above relation applies only 
in the flowing states, i.e. for values of the shear stress 
greater than the yield stress a c . The constant k is the 
Bagnold parameter introduced earlier in Eq. 1, while 
the second constant (3, having dimensions of a viscosity 
measures the relative importance of the long-lived con- 
tacts to the transient collisions. Clearly, this modified 
constitutive relation leads to a new velocity profile v x (z) 
for the material below the plug. In Fig. [3] we plot two 
best fits to the velocity profile in the flowing state. We 
have chosen data corresponding to A — 0.8, 9 = 22° that 
shows a typical example of the flow profile in the strongly 
cohesive limit. The curve (+) is the fit of the data to the 
Bagnold velocity profile where we have used the least- 
squares method. Note that one cannot simultaneously 
fit this initial slope of the velocity profile and match the 
curvature of the data. The modified Bagnold relation 
gives the best fit to the data. 

Simply demonstrating a better fit with the modified 
Bagnold form is not, in itself, conclusive evidence of the 
breakdown of Bagnold scaling as one expects a better fit 
from a model with an extra adjustable parameter. To 
better justify our conjecture we first examine the dimen- 
sionless ratio of the viscous stress proportional to /3j to 
the Bagnold stress, K7 2 : 



(14) 
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FIG. 5: (color online) Velocity profile for A = 0.8, 9 = 22° 
showing the difference in the standard Bagnold fit (blue +) 
and the modified Bagnold fit (red □). 
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FIG. 7: fi vs k n for cohesionless granular flows for H — lOOd 
and 9 = 22°. Note that the Bagnold constitutive law (fi = 
0) desribes the data better as the stiffness of the particles 
increases. 
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FIG. 6: Q vs A for 9 = 22° (closed symbols), 24° (open 
symbols) at three different heights in the flowing slab (o for 
z = H/A, □ for z = H/2, and A for z = 3H/A). Q is a 
measure of the importance of the linear term in a xz , as A 
increases so does the importance of this linear term 



For the Bagnold hypothesis, f2 vanishes. In essence we 
can determine the extent of the breakdown of such Bag- 
nold scaling by extracting the ht parameters K, /?, as well 
as the shear rate from velocity profiles as shown in Fig. [5] 
in order to compute the value of f2 from our numerical 
data. This ratio can be computed for a variety of samples 
having different tilt angle and different cohesive energies. 
In addition the ratio can be computed at different heights 
within the flowing layer of a given sample. 

In Fig. El we show a plot of fl vs. A. Examination of 
this figure reveals three trends. First, we note that for 
any height z within the flowing part of the slab and for 
any tilt angle 6, fl increases with increasing cohesive en- 
ergy A. The deviation for a purely Bagnold constitutive 
law as measured by Q increases with height in the slab 
for all measured angles and cohesive energies. For a given 
height in the slab and a given cohesive energy one finds 
that f2 decreases within increasing tilt angle 9. 

The increase of f2 with A for all heights and angles 
strongly suggests that the internal cohesion is primar- 
ily responsible for the breakdown of the Bagnold scal- 
ing. That the magnitude of the discrepancy between the 
observed flows and those derived from the Bagnold hy- 
pothesis depends on the height within the slab (17 vs. z at 
fixed A and 6) is consistent with our hypothesis that long- 
lived contacts cause this discrepancy. Clearly as one ap- 



proaches the plug the intergrain contacts have extremely 
long life-times and the constitutive relation for the ma- 
terial most greatly deviates from the Bagnold form. It 
is reasonable to suppose that in the region directly be- 
low the plug where the gravitational shear stress is only 
slightly greater than the yield stress of the material that 
there would be the greatest density of long-lived con- 
tacts in the flowing state. Thus, one would expect the 
greatest deviation from Bagnold scaling near the plug. 
The data clearly support this; for all tilt angles and co- 
hesive energies SI increases with height in the slab. Fi- 
nally, as the tilt angle is increased at constant cohesive 
energy, the total kinetic energy of the system is increased. 
The grains, having higher typical velocities will now have 
fewer long-lived contacts and transient forces associated 
with brief intergrain collisions will become the more dom- 
inant momentum-transfer process in the material. Thus 
the effective constituent law will appear to be closer to 
the Bagnold form. 

It is apparent from Fig. that even in the limit of no 
cohesion there remains a significant deviation from the 
original Bagnold constitutive law. This residual discrep- 
ancy can be well accounted for by slight interpenetrat- 
ibility of the particles in the simulation. The Bagnold 
constitutive relation can hold exactly only in the limit 
of hard sphere particles 30] . To test this point we have 
examined cohesionless particles of varying stiffness k n for 
H = lOOd. The data in Fig.[7]shows that the residual de- 
viations from Bagnold behavior of our cohesionless model 
granular material can be attributed to the finite stiffness 
of the constituent particles. 

To more directly test that the breakdown of the Bag- 
nold constitutive relation is due to the growth of the 
number of long-lived contacts in the flowing states with 
increasing cohesive energy we compiled a contact time 
histogram. We measured the times over which two par- 
ticles remained in contact allowing for both rolling and 
slipping at the contact. Collecting this data for a few 
representative runs having either no cohesion (A = 0) or 
strong cohesion (A = 0.8) we produced two contact time 
histograms as shown in Fig. |SJ Both histograms reflect 
the contact times in the horizontal slice of the slab be- 
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FIG. 8: Contact time histogram within the flowing part of 
the slab for cohesive energies A — 0.0 (black) and A = 0.8 
(blue) with tilt angle 9 = 22°. The shortest contact time bin 
has been removed for clarity. 



tween 70 < z < 100 over a period of lOr. This vertical 
height was chosen so that it is entirely contained within 
the flowing regime of the cohesive slab. In fact, from our 
discussion of Fig. we would expect the largest popula- 
tion of long-lived contacts to be found in the upper part of 
the flowing region near the plug. For the cohesive system 
studied using the contact time histogram, we then expect 
the greatest density of long-time contacts to occur near 
z = 200. In both samples the most common contact time 
is rather short: < ^contact — O.Olr. This first bin has 
been removed for clarity. Comparing the two histograms 
we note that although there are long-lived contacts in 
both the cohesionless and the cohesive granular material, 
the number of such contacts is dramatically increased in 
the cohesive case. As discussed above the presence of 
long-lived intergrain contacts in the cohesionless case is 
at least possibly due to the nonzero compressibility of the 
particles. The enhancement of the number of such long- 
lived contacts, particularly those having lifetimes greater 
than ~ 0.3t, with cohesion qualitatively supports our hy- 
pothesis that such contacts are primarily responsible for 
the breakdown of Bagnold scaling in the cohesive gran- 
ular media. The relationship between this fraction of 
long-lived contacts and a quantitative model for the pa- 
rameter will be explored in a future publication 31]. 



IV. SUMMARY 

We have studied the effects of cohesion in steady-state 
3D chute flow. Cohesive granular materials generically 
form two different flow regimes: a solid region (the plug) 
extending below the free surface characterized by a van- 
ishing shear rate of strain and a flowing region below the 
plug. The width of this plug region is independent on 
total depth of the slab but dependent on the tilt angle, 
cohesive energy, the mass density of the granular mate- 
rial and is accounted for by postulating the appearance 



of a finite yield stress in the material. The value of this 
yield stress can be reasonably determined from estimates 
of the mean cohesive force within the slab and the mean 
number of such cohesive contacts. 

We find that in the flowing state below the plug the 
velocity profile is similar to the Bagnold velocity pro- 
file: v ex z 3 / 2 . Nevertheless there remain significant dis- 
crepancies between the observed flow profile and the pre- 
dicted Bagnold form. In light of these deviations from 
the Bagnold form, we suggest a modified form of the Bag- 
nold constitutive law that combines two essentially inde- 
pendent methods of momentum transfer in the flowing 
state. This modified Bagnold relation accounts for the 
usual collisiona! momentum transfer between grains as 
well as for momentum transfer due to the cohesive forces 
acting through the long-term contacts between grains in 
the flowing state. This modified profile fits the data more 
closely. Extracting from these fits to the data the ratio of 
the stress transferred via each method, we find that this 
ratio exhibits a few reasonable trends. The effect of the 
long-term contacts grow with cohesive energy and within 
a single slab as one approaches the plug. This ratio de- 
creases with larger tilt angles; faster flowing slabs should 
more efficiently break-up any long-lived particle clusters. 

To further test this hypothesis we directly measured 
a contact time histogram within the flowing part of the 
slab. A comparison between noncohesive and cohesive 
systems indeed shows a larger number of long-time con- 
tacts in the cohesive system. This result qualitatively 
supports our hypothesis that these contacts are respon- 
sible for the breakdown of the Bagnold constitutive law. 
To provide a true quantitative model for the modified 
Bagnold constitutive law that we propose it is necessary 
to calculate from a more microscopic model the stresses 
transmitted by these long-lived contacts. In order to do 
so we must determine the spatial correlations between 
such contacts. These contacts may, in fact, form spa- 
tially correlated clusters or chains that span stream-lines 
in the flow or perhaps represent randomly distributed 
pairs of particles (dimers) that remain bound for meso- 
scopic periods during the flow. Clearly, at a fixed density 
of such long-lived contacts the stress transmission due to 
these contacts in the former case would be significantly 
larger than in the latter [3l| . 

Understanding the constitutive relation of cohesive 
granular materials is clearly of fundamental importance 
to the study of granular flows in a geophysical context as 
well as to the handling of granular materials in industry 
where small amounts of a wetting fluid create adhesion 
contacts between the grains. Interestingly, this work sug- 
gests that studying granular flows in the presence of small 
amounts of cohesion allows one to break the expected 
Bagnold scaling in a controllable and computationally 
efficient manner. It is clear from this work and others 
[Til l30| that the non-zero compliance of the grains leads 
to measurable deviations from the Bagnold scaling; the 
Bagnold law holds in the limit of hard spheres. Studying 
significantly less compliant grains, however, is computa- 
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tional difficult. The further study of cohesive granular 
materials both analytically and computationally should 
enable the exploration of granular constitutive laws for 
physically accessible systems that nearly, but not exactly 
obey the Bagnold constitutive relation. 



APPENDIX A: DERIVATION OF THE 
VELOCITY FROM THE BAGNOLD 
CONSTITUTIVE LAW 



therefore 



Oxz{z) = pgsin9(H - z), 



(A3) 



where H is the total height of the slab. Solving for v x (z) 
gives the velocity profile shown in Eg. H2l above. 

Modified Bagnold scaling is derived from the same con- 
stitutive relation used in deriving traditional Bagnold 
scaling [Eq. lAlj with the addition of a linear term in 
7 



Bagnold scaling is derived from a constitutive relation 
between the shear stress and the strain rate 



cr X z = «7 



(Al) 



where 7 = dV g^ ■ The steady-state Cauchy equation for 



o xz is 



da x 



pg sin6>, 



«7 



(3j. 



(A4) 



Solving this with the same Cauchy equation from 
Eq. (|A3fl results in 



v x (z) 



(G 2 + cH) 3/2 - (G 2 + c(H - z)) 



3/2' 



( A2) where c = pg sin 9 and G 



-Gz, 
(A5) 



O. Pouliquen and F. Chevoir, Physique 3, 163 (2002). 
A. Suzuki and T. Tanaka, Ind. Eng. Chem. Fundam. 10, 
84 (1971). 

D. A. Augenstein and R. Hogg, Powd. Tech. 19, 205 
(1978). 

T. J. Drake, J. Geophys. Res. 95, 8681 (1990). 

C. Ancey, P. Coussot, and P. Evesque, Mech. Cones- Fric. 

Matter 1, 385 (1996). 

O. Pouliquen, Phys. Fluids 11, 542 (1999). 

E. Azanza, F. Chevoir, and P. Moucheront, J. Fluid 
Mech. 400, 199 (1999). 

C. Ancey, Phys. Rev. E 65, 011304 (2002). 

O. Pouliquen, Phys. Rev. Lett. 93, 248001 (2004). 
O.R. Walton, Mech. Mater. 16, 239 (1993). 
X.M. Zheng and J.M. Hill, Powd. Tech. 86, 219 (1996). 
S. Dippel and D.E. Wolf, Comp. Phys. Comm. 121, 284 

(1999) . 

D. M. Hanes and O.R. Walton, Powder Technol. 109, 133 

(2000) . 

L.E. Silbert, D. Erta§, G.S. Crest, T.C. Halsey, D. 
Levine, and S.J. Plimpton, Phys. Rev. E 64, 051302 

(2001) . 

L.E. Silbert, G.S. Crest, S.J. Plimpton, and D. Levine, 
Phys. Fluids 14, 2637 (2002). 

L.E. Silbert, J.W. Landry, and G.S. Crest, Phys. Fluids 
15, 1 (2003). 

P.K. Haff, J. Fluid Mech. 134, 401 (1983). 



[18] S.B. Savage, Advances in Micromechanics of Granular 
Materials eds. J.T. Jenkins and M. Satake (Elsevier, New 
York) (1992). 

[19] S.B. Savage, J. Fluid Mech. 377, 1 (1998). 

[20] J.T. Jenkins and S.B. Savage J. Fluid Mech. 130, 187 
(1983). 

[21] W. Losert, L. Bocquet, T.C. Lubensky, and J.P. Gollub, 

Phys. Rev. Lett. 85, 1428 (2000). 
[22] L. Bocquet, W. Losert, D. Schalk, T.C. Lubensky, and 

J.P. Gollub, Phys. Rev. E 65, 011307 (2001). 
[23] L. Bocquet, J. Errami, and T.C. Lubensky, Phys. Rev. 

Lett. 89, 184301 (2002). 
[24] P. Mills, D. Loggia, and M. Tixier, Eur. Phys. J. E 1, 5 

(2000). 

[25] M.Y. Louge, Phys. Rev. E 67, 061303 (2003). 
[26] M.Z. Bazant, unpublished (2005). 

[27] R.A. Bagnold, Proc. R. Soc. London, Ser A 225, 49 

(1954); 249, 235 (1956). 
[28] P.A. Cundall and O.D.L. Strack, Geotechnique 29, 47 

(1979). 

[29] D. Erta§ and T.C. Halsey, Europ hys. Lett. 60, 931 
(2002); T.C. Halsey and D. Erta§, |cond-ma t/0506170 
(2005). 

[30] G . Lois, A . Lemaitre, and J.M. Carlson, 

cond-mat/0501535 (2005). 
[31] R. Brewster and A.J. Levine, in preparation (2005). 



